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ABSTRACT 

We present an estimate of the third integral of motion for axisymmetric three-dimensional potentials. This estimate is based on a 
Stackel approximation and is explicitly written as a function of the potential. We tested this scheme for the Besanfon Galactic model 
and two other disc-halo models and find that orbits of disc stars have an accurately conserved third quasi integral. The accuracy ranges 
from of 0.1% to 1% for heights varying from z = 0 kpc to z = 6kpc and Galactocentric radii R from 5 to 15kpc. We also tested the 
usefulness of this quasi integral in analytic distribution functions of disc stellar populations: we show that the distribution function 
remains approximately stationary and that it allows to recover the potential and forces by applying Jeans equations to its moments. 
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1. Introduction 

Many approaches exist for the galactic modelling of stellar dy¬ 
namics. With the availability of a large amount of new accurate 
data for the velocity and position of stars in the extended solar 
neighbourhood, the development of the most precise tools for 
analysing the kinematics of stars is becoming inescapable. 

Numerical modelling techniques are the most direct ones, 
such as for the Scharwschild method, the made-to-measure 
method ( |Syer & Tremaine|1996[ ), or of course ab initio N-body 
and hydr odynamical simulations ([Renaud et al.|2013|l. Spectral 
analysis ( |Papaphilippou & Laskar|1998[|Valluri & Merritt| 1998) 1 
also offers useful numerical tools for recognising and identifying 
resonances. Numerical resolution of the collisionless Boltzmann 
equation is another direct way to model, but it is still limited by 
available numerical resources (|Yoshikawa et al.|2013l|Colombi| 
et al.|2015F 

Analytical techniques are more explicit but imply analytical 
approximations. In the case of axisymmetric galactic potentials, 
it is known that orbits are generally constrained by an effective 
third integral of motion in addition to the energy and angular mo¬ 
mentum. From the Jeans theorem, therefore, equilibrium models 
of axisymmetric galaxies should be represented by distribution 
functions that only depend on these three integrals. The mod¬ 
elling of stellar distribution functions can thus be achieved if the 
effective third integral can be approximated analytically or nu¬ 
merically. 

A long list of works in this direction already exists, and 
very useful bibliographies may b e found in|de Zeeuw ( 1985|l, 
de Zeeuw & Lynden-Bell (|1985 l,|De Bruyne et al. ( |2000 l, Bin-| 


fitting the true potential with a Stackel potential. This allows 
nearly exact modelling of orbits and integrals of motion to be 
obtained in the immediate neighbourhood of the position con¬ 
sidered. 

Global modelling, either of the potential or of the orbits over 
a predefined phase-space volume, has opened the way to many 
different technicalities over the past 50 years. In each case the 
goal was to use a Stackel third integral as an approximate of 
the third quasi integral of the fitted potential, as in Wayman 


( | 1 959| l, [van de Hufst (|1962[), Ollongren \ \962\, Hori (| 1962^, and 


more recently, Manabe (|1979[ ), Batsleer & Dejonghe ( 1994| l, or 
Famaey & Dejonghe|p()03lk 


ney ( 2012[ ), and Sanders & Binney| ( |2014| l, among others. We 
concentrate here on the use of Stackel potentials that have been 
shown in many cases to be efficient for modelling different fam¬ 
ilies of orbits. Published works on this subject present a wide 
variety of approaches, and we can distinguish between the local 
and global approaches. The local modelling consists in locally 


A recent novel application of the Stackel potential approxi¬ 
mations has been to model the local stellar kinematics and the 
gravitational potential in the solar neighbourhood at large z dis¬ 
tances up to 1-2 kpc from the galactic plane ( |Bienayme et al.| 
|2014[[Prffl et al.|2014| l. These recent studies put strict constraints 
on the vertical variation in the gravitational potential and in the 
local density of the dark matter halo. Recently, [Sanders & Bin-| 
jneyj ( |2014| l have proposed a global Stackel fitting of 3D poten¬ 
tials, and they give algebraic expressions that allow recovery of 
the integrals of motion, expressed as action variables, and they 
present numerical applications. 

In the present paper, we proceed to a Stackel potential fitting 
by using a simple expression for the integral of motion that ex¬ 
plicitly depends on the potential. Our study has similarities with 
the works recently published by|Binne^(|2012|) or byjSanders &| 


Binney 


([2014 1, proposing different formulations of a third inte¬ 
gral. Although there are also advantages in working with action 
integrals, the present approach is, however, much simpler and 
more straightforward to apply. 

The paper proceeds as follows. In Section 2, we give a new 
expression for the third integral, and its derivation is described in 
the Appendix. In Section 3 we examine, for three potentials, the 
constancy of this third integral along orbits. Section 4 shows, in 
the case of the Besan 9 on Galactic model, how moments of dis- 
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tribution functions based on this third integral are in accordance 
with Jeans equations. Section 5 finally considers the application 
to the collisionless Boltzmann equation. 


2. Quasi integral of motion 

We let V{R,z) be an axisymmetric gravitational potential. Be¬ 
sides the energy E and the vertical component of the angular 
momentum L-, which are integrals of motion, we dehne here, as 
an approximate third integral of motion. 




\L^-Li 



( 1 ) 


where L and L, are the total and vertical angular momenta, zo 
a fixed parameter at fixed E and L^, the vertical velocity, and 
with 


>?(/?, z) = -[y(7;,z)-l/(Vd,0)] 



( 2 ) 


3. Testing the conservation of the integral along 
orbits 

We examined the conservation of the quasi-integral of motion 
along orbits of stars corresponding to disc stellar populations. 
For these disc components, stellar motions have restricted oscil¬ 
lations in the radial and vertical directions. For the smallest oscil¬ 
lations, the phase-space domain explored by stars is the closest 
to a Stackel potential, and the quasi integral tends to be constant. 

We considered stars with identical angular momentum 
(with Ec the energy of the corresponding circular orbit) and iden¬ 
tical energy E (we dehne AE - E - Ec). For each star, we com¬ 
puted the mean value of Is and its dispersion ct/j from ~10 000 
steps along the orbit, and we determined the value of zo that 
minimizes ajs- To ensure that cr^ is determined sufficiently well, 
the total time integration for each orbit is about 300 dynamical 
times (i.e. 300 galactic rotations). The numerical integration is 
a Runge-Kutta-Fehlberg of order 7(8) ( Fehlberg |1968| ). We hnd 
that stellar orbits with the same values of L~ and AE all have 
<T[s minimum for similar values of zo- For different energies and 
angular momenta, the adjusted zo will be different. 

To allow an easier presentation of our results, we normalized 
L as 


and 

{(R^ + - Zo) + 2 V + 4 z2 . 


(3) 






(4) 


This approximate integral /j is an exact integral in the case 
of Stackel potentials expressed in an ellipsoidal coordinate sys¬ 
tem of focus Zo, which is thus a known quantity. Its derivation is 
explained in the Appendix. 

The only free quantity for determining /,, when V{R,z) is 
known, is the parameter zo- Here, it will be numerically adjusted 
by minimizing ctjs the dispersion of the quasi integral Is along all 
the orbits with the same energy and vertical angular momentum. 
Thus Zo will itself be a function of the two integrals E and L^. 

In a preliminary analysis, which was not developed further, 
we used the generic differential equation of Stackel potentials 
that can be rewritten to determine zo locally as a function of the 
second derivatives of the potential at any positions (see Eq. 10 
Bienayme|2009[ de Zeeuw|1985|l. However, we hnd out that 


such a procedure gives us signihcantly fewer accurate results 
than htting a single zo for a given orbit or family of orbits. Here, 
we look for the zo value that minimizes the dispersion of Is along 
each orbit of a family of orbits with the same E and L^. We note, 
however, that for any (£’,L,)-family of orbits and with Zsheii being 
the maximum z-extent of the shell orbit, the value of zo that gives 
the best ht is close to the value of zo obtained by applying Eq. 
10 of Bienayme ( 2009 | l at position (Rc, Zsheiil^) (with Rc(Lc) the 
radius of the circular orbit with angular momentum L,). 

The degree of the approximation of the quasi integral Z, can 
be controlled in several ways: ( 1 ) inspection of the surfaces of 
section, ( 2 ) conservation of orbital weights and of the spatial 
density, and (3) conservation of Is along the orbits to validate 
the labelling of orbits. The conservation of orbital weights and 
the conservation of Is are the criteria that can be best expressed 
numerically. 

In the following sections, we test the degree of conser¬ 
vation of this quasi integral along orbits in the cases of three 
non-Stackel potentials. We also determine its efficiency to build 
stationary stellar disc distribution functions and to measure the 
potential from the Jeans equations. 


still an integral of motion. 1^-0 corresponds to orbits conhned 
within the galactic plane, while /3 ~ 1 corresponds to shell orbits 
with minimum radial variations when they cross the plane z = 0 . 

In the following sections, we consider three potentials: a log¬ 
arithmic potential with a disc component, a flattened logarithmic 
potential, and a potential similar to the Besan 9 on Galactic model 
potential. In all the three cases, the constancy of the third inte¬ 
gral is satisfied at 0.2 to 0.5 per cent for orbits corresponding to 
thin or thick stars and to a 1 to 2 per cent for orbits with larger 
vertical oscillations. 

Table 1. Disc-halo logarithmic potential. Maximum and median val¬ 
ues of the quasi-integral dispersion cr /3 (i.e. the relative error) for orbits 
with different energies AE: best fit of zo, maximal height Zmax above the 
galactic plane for the shell orbits. 


AE 

Zo 

^max 

cr /3 maximum 

(T /3 median 

0.05 

3.3 

2.4 

0.008 

0.004 

0.1 

3.3 

3.7 

0.006 

0.003 

0.2 

3.3 

6.1 

0.021 

0.005 

0.5 

3.3 

12.9 

0.012 

0.006 


3.1. Logarithmic potentiai of Myamoto-Nagai type 

The logarithmic potential of Myamoto-Nagai type 


log + 


-I- z^j j 


(5) 


has a positive density ( |Zotos||2OTT i. It has a spherical halo and 
a disc component (with a + b the halo core radius and b the 
disc thickness). This potential presents analogy with the family 
of density-potential pairs of logarithmic Myamoto-Nagai type 
given by|Bienayme|(|200^. 
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We consider the orbits with the angular momentum L^—8.5 
(n<xi = l, a-0.5, b-0.2, Rc =8.5) and with E - Ec + lS.E. The initial 
conditions of the computed orbits are vr - 0 , with Rinuiai equally 
spaced. These initial conditions do not include resonant orbits 
that do not cross the galactic plane perpendicularly. 

Table 1 summarizes the results. For each energy AE, the ad¬ 
justed zo that minimizes cr /3 is given. The maximum vertical ex¬ 
tension Zmax corresponds to an orbit with Ij, close to 1. The dis¬ 
persions of 73 always remain small, the mean dispersion being 
smaller than one per cent. However, for some resonant orbits, 
the dispersion is larger, of the order of 2 per cent. We also note 
that at large energies corresponding to halo stars the variation of 
73 still remains small. This must be linked to the spherical shape 
of the potential at large z, easily modelled by a Stackel potential. 

In the next section we consider a different requirement using 
a flattened potential. 



R 


3.2. Logarithmic potentiai 

To examine the impact of the flattening of the dark halo on the 
dispersions cr/ 3 , we consider the traditional logarithmic potential 
|Richston^ ( |1980| l: 




( 6 ) 


We set Vc=l and q - 0.8 to correspond to a flattening of the 
isodensity qp - 0.65. 

We computed orbits with - 8.5 {Rc - 8.5, Ec - 1). Table 
2 summarizes the results. Surfaces of section for Lj = 1, AE - 
0.05, and 0.5 are shown in Fig. 9 of Bienayme & Traven ( 2013| l, 
and we note the lack of resonant orbits. This may explain the 
excellent conservation of h at a level of 0.2 per cent for stellar 
disc orbits. Orbits with larger extension up to Zma.-c-9 (7i=0.4) 
have 73 conserved at 2 per cent. 

The inspection (see Fig. [^l of the surface of section and of 
the meridional projection of three orbits (AT: =0.2 and L,=8.5) al¬ 
lows visualization of the relative agreement between the nearly 
exact numerical orbits and the contours obtained from the ap¬ 
proximated third integral. For the meridional projection, the red 
dashed lines are ellipsoids and hyperboloids from the ellipsoidal 
coordinates with zo=5.9 (lines are drawn from the corner of the 
orbit envelopes). The blue continuous lines in Figure [1] (middle 
panel) are the envelopes that are determined semi-analytically 
from E, Lj and the quasi integral 73 that are explicitly known; 
we determine surfaces of section at various fixed R (resp. fixed z) 
and find the z (resp. R) extension limits of the orbit from the con¬ 
dition dIj,ldvr-Q (resp. dIj,fdv^=Q). We note small differences 
between the envelopes of numerical orbits (top panel) and semi- 
analytical envelopes (middle panel). We also note that the analyt¬ 
ical hypersurfaces defined from E, do not have the topology 
of a torus everywhere. The bottom panel of Figure [T shows a 
section of the phase space (z= 0 ) for the same three orbits from 
the numerically computed orbits and from the analytical section 
given by Eq[T] If lic= 200 km/s, the maximum vr of these three 
orbits is ~120, 50, and 20 km/s. 



R 



R 


Fig. 1. Top panel: Meridional projection for three orbits within the 
logarithmic potential (Sect. 3.2) with 7.-=8.5 and A£=0.2. Zero veloc¬ 
ity curve and some coordinate curves of the elliptic coordinates (with 
Zo=5.9) are drawn. Middle panel: Analytical determination of the en¬ 
velopes of the orbits. Bottom panel: Surfaces of section for the same 
three orbits. Black crosses: numerically computed orbits. Red lines: the 
corresponding sections obtained from Eqs |l|3| Blue line: surface of sec¬ 
tion of the orbit confined in the mid-plane (i.e. U;;=0). 


3.3. The Besangon Gaiactic mass modei 

We apply our test to the gravitational potential of the Besan 9 on 
Galactic model (BGM) ( [Bienayme et al.|1987[|Robin et al.|2003[ 
Czekaj et al.||2014} , which is a more realistic one than the two 
previous ones, although we limited our analysis to its axisym¬ 
metric version. This model has a dynamical consistency at the 


solar Galactic radius position in the sense that the thickness of 
the stellar components are constrained by the potential and by 
the vertical velocity dispersions. However, outside of the solar 
neighbourhood, the vertical kinematics and the thickness of stel¬ 
lar components are constrained by observations and not by a dy¬ 
namical consistency of the model. The motivation of this work is 
to improve the BGM dynamical consistency by using stationary 
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Fig. 2. Left: quasi integral versus Zmax for stars with AE = 6400 and f?c(L-)=8500pc for the BGM potential. Errors bars are cr/ 3 . Right: idem 
with AE = 12800. 





■.. 


° 

0.03 

- 

0. 04 




0. 03 

• : 

0.02 


m 

0 

n 


g 

- 

g 

s 

0.02 

- 

ui 



° 0 




] 


0 \ • 

0. 01 

0 “ ft 




° e 


) 

0 


0 

- - y 


0 1000 2000 3000 0 2000 4000 6000 

Zmax (pc) Zmax (pc) 


Fig. 3. Left: Dispersion <T /3 of the integral along orbits versus Zmax for stars with AE = 6400 and Rc(L,)=8500 for the BGM potential. Right: 
the same with AE = 12800. 


Table 2. Logarithmic potential (see legend of Table 1). 


AE 

zo 

Zmax 

cr /3 maximum 

CT/a median 

0.04 

5.4 

2.0 

0.004 

0.002 

0.1 

5.5 

3.3 

0.011 

0.005 

0.2 

5.9 

5.2 

0.020 

0.010 

0.4 

7.0 

8.8 

0.040 

0.019 

0.8 

10 . 

17.5 

0.08 

0.04 


distribution functions to model the stellar kinematics. Besides 
the integrals E and L,, we therefore require an approximate third 
integral to describe the stellar kinematics. 

We determined the gravitational potential of the Besan 9 on 
Galactic model according to the characteristics of its compo¬ 
nents described in [Robin et al.| ( |200^ and jCzekaj et al.| ( |2014| l. 
The dark matter component is represented by a spherical com¬ 
ponent. The orbits with large vertical extensions have shapes 
that are mainly determined by the spherical halo. Since the exact 
shape of the dark halo remains uncertain in practice and is a con¬ 
troversial question, we consider a less favourable situation using 
a flattened dark matter component. 

Article number, page 4 of|^ 


The dark matter halo density used in the BGM is 


Pdm ~ 


1 

1 + (R/Rcore)^ ■ 


(7) 


We flatten the halo, replacing the quantity R with in 

the expression of the potential. With Rcore=2700pc and q - 0.8, 
the flattening of the density is qp - 0.62. (For this value of q, 
negative densities only occur in a region far from our domain of 
interest.) 

We also modify the point-mass bulge ( [Bienayme et al.|1987| l 
using the potential law 0)^ = with Rhuibe = 

2000 pc. 

We summarize the results for orbits with corresponding to 
Rc =8500 pc in Table 3. Results are listed by families of different 
energies AE. Figure]^ shows /a versus Zmax, the maximum verti¬ 
cal extension of orbits in cases A£’=6400 or 12800 and Rc =8500 
pc. Error bars are the dispersion of along each orbit. The dis¬ 
persion of /a is zero for orbits confined in the mid-plane. Shell 
orbits have /a maximum and cr/a close to a few 10“"^ (Figure 
|3 I. This is unexpected since there is no reason that the potential 
along the ’thin’ orbits should have exactly the Stackel form. On 
the other hand, when Zmax~^000pc, the dispersion of /a is max¬ 
imum and increases sharply for all energies (see Figure [^, also 
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corresponding in Figure|^to a change of shape of the distribution 
of points. This is related to the presence of (1,1) resonant orbits 
at this vertical height that are poorly modelled by our Stackel ad¬ 
justment. For these orbits the dispersion of the quasi integral is 
less than 5 per cent. The median dispersion of for orbits with 
Zmax smaller than 2 kpc remains very small ~ 0 . 002 . 

Table 3. Besanfon Galactic model with ijr=0.8 and l?c=8.5 kpc (see leg¬ 
end of Table 1), where is the maximum vertical velocity. 


AE 

^max 

Zo 

Zmax 

cr /3 max. 

cr /3 med. 

km^.s~^ 

km.s“^ 

kpc 

kpc 



400 

28.3 

4.8 

0.42 

0.001 

0.0004 

1600 

56.6 

4.4 

1.0 

0.002 

0.0012 

3200 

80 

4.3 

1.7 

0.004 

0.0025 

6400 

113 

4.5 

3.5 

0.047 

0.0043 

12800 

160 

5.1 

6.1 

0.032 

0.0070 


We performed the Stackel adjustment for other galactic radii 
Rc from 1.5 to 15 kpc. For each pair of values {E - Ec + AE, L^), 
we determine the zq value that minimizes (T/ 3 . Figure j^plots the 
fitted zo (crosses) for A£’=100,200,400...12800. Lines of con¬ 
stant AE are plotted. We obtain a tabulation of zo that gives us a 
function of two integrals of motion zoiE, L^) that can be used to 
build a distribution function for stellar populations. 

The median dispersions of remain low except at galactic 
radii smaller than 3 kc where many resonant orbits dominate the 
phase space. Figure|^plots the histograms of cr /3 for R,. from 3.5 
to 15.5 kpc. 



Rcirc (pc) 


Fig. 4. BGM potential: zo best fit versus Rcirc(Lz) and AE. Black lines 
A£=100 and 200; red lines: 400 and 800; green lines: 1600 and 3200; 
blue lines: 64000 and 128000. 


4. Distribution function and Jeans equations 


We tested the efficiency of using the quasi integral to model 
the distribution function of disc stars using a Shu distribution 
function (DF) generalized to 3D axisymmetric potentials. It 
writes as (here, we correct a typographic error in IBienaymel 

|T^ : 


f(E,L„h) 


2D.(R,) 2(L,) 
iTTKiRc) O-l 


(E - Ecirc) 


X 


1 1 


exp 

- (£ - Edrc) 

' 1 1 ' 
2 2 

h 



' ^R . 




Fig. 5. Histograms of E dispersion along orbits for Zmax intervals de¬ 
limited by 0 , 1000 pc, 2000 pc, and beyond (respectively black, dotted 
red, dashed blue lines). For each of the three Zmax intervals, the medians 
of 0-/3 are respectively 5. 10“'*, 1.5 10“^, and 5. 10“^. 


with Rc(E^) the radius of the circular orbit with the angu¬ 
lar momentum L^, Q. the angular velocity, k the epicyclic 
frequency, and Edrc the energy of a circular orbiting star 
at radius Rc- For sufficiently small velocity dispersions, the 
number density distribution, £(Lj.) = Sq exp(-^c/^v), is close 
to 2(7?) = So exp(-R/Ry). We set constant the parameters 
and /?v=2.5kpc, which is close to the scale length 
of the number density distribution. The DF allows us to re¬ 
produce the triaxiality and tilt of the velocity ellipsoid and to 
model nearly exponential density disc distribution. It could 
be easily modified to reproduce any reasonable radial density 
law. We set a constant velocity dispersion {Rg.^_-co) for a 
thin disc (ctr, crj)=(40km/s,20km/s) and for a thick disc 
(ctr, crj)=(60km/s,40km/s). 


The distribution function can also be written in a more read¬ 
able form as 


f(E,L.,h) = g(L-)e.xp 





2 

exp 

2 



O^zJ 


(9) 


where £r and £- are integrals of motion that can be easily de¬ 
duced from Eq. 


£r =(£-£„>,)(!- 73 ) 

a HE-Edrc) h- 


They are respectively related to the amount of radial or vertical 
motions that can be controlled versus the parameters ctr and 
cTj. Within a Stackel potential and for an orbit with angular 
momentum they are respectively at position {RdL^), z = 0 ) 
the radial and vertical kinetic energies. Shell orbits have £r=0. 


We determine the moments of the distribution function (Eq. 
1^ and recover the radial and vertical forces, Kr and K^, from 
the Jeans equations for a stationary axisymmetric potential. The 
moments are linked through the Jeans equations where the time 
derivative are not exactly zero since 73 is just an approximate 
integral; 


dvvR d — d _^ 


^2 2 

R 




= 0 , 


( 11 ) 
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Fig. 6. (Left) Dark lines: radial forces Kj; of the BGM at several z above the galactic plane (0.5, 1, 2, and 3 kpc bottom to top). Thin green lines: 
Kn recovered from a Jeans equation in the case of a thick disc. Thin red lines: the case of a thin disc. (Right) Relative errors on the recovered 
radial force Kj; versus R at various z (0.5, 1, 2, and 3 kpc) (resp. green, black, red, and blue lines) for thin disc DF (dotted lines) and thick disc DF 
(continuous lines). 


dvV^ d _ ^ / “Tn / r. /1-.X 

-aT * * s'"'' * ■' T" & r “■ 


By considering that is close to an exact integral and neglecting 
the time derivative terms, we transform the Jeans equations as 


d -r d _ 




R 




= 0 , 


(13) 


d _ ^ 1 


= 0 . 


(14) 



z (pc) 


Thus, under the assumption that the approximate integral 
is exact, Eqs. T3p4 are used to obtain an estimate of the radial 


and vertical forces, Kr 


\aR)e 


K, 


•(^) ^. They can be 


compared to the exact forces to test the efficiency to recover the 
galactic potential using our Stackel integral in a non-Stackel po¬ 
tential. Figure|^shows the radial forces Kr of the BGM potential 
and its relative error at several z when recovered from Eqs [T3p4| 

Within the domain R = 3 to 16kc and z smaller than 3 kpc, 
the relative error on the radial force Kr is smaller than one per 
cent for R > 6 kpc. The relative error is ten per cent at /? = 3 kpc 
and increases below R = 3 kpc. 

Figure [7] shows the exact and recovered K^ force versus z at 
several R (3.5, 4.5, 6.5, 8.5, 12.5). At large R > 6 kpc, the K^ 
forces are accurately recovered for the thick disc up to 6 kpc and 
up to 3 kpc for the thin disc. The vertical force K^ at the solar 
position J?=8kpc is remarkably well recovered at 0.5% up to 
6 kpc for the thick disc (7 scale heights) (Figure]^. For the thin 
disc, the accuracy is lost above 3 kpc (corresponding however to 
ten scale heights of the thin disc). At lower radius R from 3.5 to 
5, the K^ force is only approximately recovered up to z-2 kpc. 

For comparison (Figure]^, we also plot the K^ force recov¬ 
ered neglecting the cross term < vrv^ > in Jeans equations, a 
classical assumption valid at low z. We note that this hypothe¬ 
sis remains valid up to 500 pc for the thin disc. At higher z, the 
recovered force diverges quickly from the exact one. 


Fig. 7. Vertical force K^ versus z at R=(3.5,4.5,6.5,8.5,12.5 kpc) (top to 
bottom) from the BGM (black lines) and recovered K. for a thick disc 
DF (blue dotted lines). 



z (pci 


Fig. 8. Top : Vertical force versus z at R=8500 pc (black line) and 
recovered from simplified Jeans equation assuming the separability 
of vertical and radial motions: red dotted line for a thin disc DF and blue 
dashed lines for a thick disc DF. 
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5. Collisionless Boltzmann equation 

The stationarity of the Jeans equation is a necessary condition 
for validating the degree of stationarity of a distribution func¬ 
tion, but this is not a sufficient condition. We should also have to 
test the stationarity of all the other moments of the collisionless 
Boltzmann equation (CBE), since it is easy to find solutions for 
the Jeans equations that are not solutions of the CBE, so the 
Jeans equation test can be more optimistic than the indications 
solely obtained from the conservation of the integral of motion. 
Eor this reason we examine hereafter the stationarity of the CBE. 


5.1. First estimate 


We estimate the stationarity of the distribution function built 
with the approximated third integral by looking at the time vari¬ 
ation of the DE (Eq. over a dynamical time (~ an orbit revo¬ 
lution). We have 


din/ d\nf 


dt 


df 


+ [ln/,i/] =0 


or 


51n/ ^E 
dt crl 




dt 


(15) 


(16) 


and the relative variation of / over a dynamical or longer time 
is determined by the variation 0-/3 of the quasi integral / along 
orbits: 


/ 


AE 


cr^ 


■0-/3. 


(17) 


Orbits with small radial or vertical amplitudes (AE and cr /3 
small) have the smallest variations in /. Thus the thinnest discs 
are the most accurately modelled. 

The stationarity decreases quickly with z since we have 
AEfcr^ oc It also decreases because of the significant pres¬ 
ence of resonant orbits between 1 to 2 kpc. However, at the solar 
Galactic radius Rq, orbits within the BGM with large vertical 
amplitude remain correctly modelled with a Stackel potential 
thanks to small cr /3 of the order of 0.005. In this situation, the 
distribution function remains stationary at many scale heights 
above the galactic plane. Eor the thick disc, the stationarity (Eq. 
171 is about one per cent at six scale heights (z ~ 6 kpc). Eor the 
thin disc, it becomes insufficient at about ten scale heights ( z ~ 
3 kpc). 

The time derivative in Eq[^is the time variation from an ini¬ 
tial position of stars in phase space at f =0 given by the initial con¬ 
dition of Eq. How divergent are the orbits evolving from this 
supposedly near-equilibrium initial condition? We know that for 
potentials similar to the three potentials considered in this paper, 
the orbits are essentially regular and an effective third integral 
must exist, which could be eventually estimated with a higher 
accuracy using, for instance, high -order polynomials (|Bienayme| 
|& Traven|2013| l or a torus fitting ( [Sanders & Binney|2014| l. This 
would, however, not be true in cases of significant presence of 
ergodicity, for instance close to the corotation of a barred rotat¬ 
ing potential. It turns out that our approximate integrals oscillate 
along each orbit around a mean value, that a typical periodicity 
of these oscillations is the dynamical time, and that the ampli¬ 
tude of these oscillations is of the order of cr/ 3 . This is illustrated 
in Eig.|^ where the maxima of d/ 3 /df are of the order of cri^ltciyn 
with tdyn ~ 6 . 



Fig. 9. d/ 3 /dt: time derivative of the quasi integral within the logarith¬ 
mic potential for the three orbits shown in Fig[2 (black h = 0.16, red 
/ =0.88, blue /3 = 1.06). 


5.2. Second estimate 

It is also useful to consider the non-stationarity as a shift in po¬ 
sitions or velocities of the modelled DE relative to a stationary 
DE. A simple ID dimensional analysis allows it to be illustrated. 
With a quadratic potential O = az^ and a stationary DE as 

fstat. = exp (-(/(z) -H vl/2)/cr^) , (18) 


the shifted DE in velocity by a factor vq is 


/ = exp 


1 



(vz - I 


(19) 


Such a DE could be considered satisfying if the shift is small 
even if we notice that the relative errors on / increase in the tails 
of the distribution. 


We deduce from the CBE, 


d\nf 

dt 



( 20 ) 


which combined with Eq. gives the variation in / over a dy¬ 
namical time tdyn- 


0-2 0-2 tdyn ' 


( 21 ) 


With a quadratic potential and h, the scale height of the disc 
population, we estimate the velocity shift: 


AE cTn z 

- - ~ —O’CTj'X 

\K,\ tdyn h 


( 22 ) 


Eor the thin disc DE at one scale height (h ~ 250 pc), the shift 
Vq is 0.2 km/s and at 4 kpc (16 scale heights), it is 10 km/s (to 
be compared to the 20 km/s of the vertical velocity dispersions) 
where the application of the Jeans equation shows that the DE is 
not stationary. At one scale height of ~ 1 kpc for the thick disc, 
the shift is Ho ~ 1 km/s. At 6 kpc, it is 7 km/s (to be compared 
to the 40 km/s of the vertical velocity dispersion), and the Jeans 
equation still validates the stationarity. 
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6. Conclusion 

In this paper we have shown that Stackel potentials can be used 
to fit a wide variety of disc stellar orbits within different axisym- 
metric potentials of disc galaxies. We also proposed a new and 
simple formulation for the third integral of motion of Galactic 
potentials, explicitly depending on the potential, which is an in¬ 
tegral known to be exact in the case of Stackel potentials. 

The quality of the fit of orbits was quantitatively measured 
by looking at the conservation of the approximate third integral 
of motion besides the energy and angular momentum. By using 
the Besan 9 on Galactic model, we showed that the third integral 
is conserved to a few thousandths in a wide volume around the 
solar neigbourhood. It is conserved to better than one per cent 
at 6 kpc above the Galactic plane at the solar position. However, 
the fit fails at low galactic radius (R smaller than 4 kpc) owing to 
the presence of a large number of resonant orbits. 

In light of the future Gaia data, we also used distribution 
functions of disc stars, depending on the three integrals, and 
checked the ability of Jeans equation to recover the gravitational 
potential. We also considered the stationarity of the CBE. 

In conclusion, Stackel potentials can be good local approx¬ 
imations of general realistic potentials. They have been used 
many times to fit the global potential of galaxies (see for instance 
[Dejonghe & de Zeeuw|1988| ) or the potential of our own Galaxy 
( Famaey & Dejonghe|2003| l. The fit of stellar orbits are also fre¬ 
quently used using Stackel potentials with local or global fitting 
( [Kent & de Zeeuw|1991| l. We showed here that they be can used 
very efficiently to build distribution function of disc stellar pop¬ 
ulations. A straightforward application will consist in using such 
modelling of DFs and in extending the dynamical consistency of 
the Besan 9 on Galactic model to describe the kinematics of disc 
stellar populations. 
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Appendix A: Quasi integral of motion 


We refer the reader to de Zeeuvt^ ( |1985| l and |de Zeeuw & Fynden 
|Bell| ( |19^ for notations and for a detailed description of prop 
erties of Stackel potentials. 

An axisymmetric Stackel potential is fully defined by two 
free functions h{A) and h(v). In the case of prolate spheroidal co¬ 
ordinates, +zo (Zq = y - ff) are the foci of the confocal ellipsoids 
and hyperboloids used to define a system of coordinates (A, v) in 
which Stackel potentials are more easily tractable. 

Fet V(R,z) - V(A,v) be the true potential that we approx¬ 
imate locally at the position (J?i,zi) = (/li,vi) with a Stackel 
potential <I)(/l, v). We assume that zo is already known. 

By definition, we have for a Stackel potential: 


d>(T, v) = - 


h(A) - h(v) 
A - V 


(A.l) 


and it is always possible to find a Sackel potential that coincides 
with any potential on the chosen coordinate surfaces A = Ai and 
V = vi. It writes as 


<l)(T,v) = 


V(A,vi)(A - vi) - V(Ai,vi)(Ai - vi) -F V(Ai,v)(Ai - v) 
A - V 


(A.2) 


Thus 

/i(A) =y(T, viXA - Vi) - y(/li, Vi)(di - Vi) + c 

h(v) - - V(Ai,v)(Ai - v) + C (A.3) 

and C an arbitrary constant. 


We can express the third integral X associated to O as 


2zl 2 

or 


1 2 


1 -F 


R 


a\ 


Zq 


2 , RzVrV^ 

'’z + — 15 — 


^0 


with 


T'(T, v) = 


(v -F y) h{A) - (A + y) h(v) 
(y -a) (A- v) 


(A.4) 


(A.5) 


(A.6) 


We fix the remaining free constant C by setting h(v = -y) = 0, 
so T* is null at z = 0 in the plane of symmetry of the potential. 
Then, evaluated at (Ti, vi), this function simplifies as 




(y(di,vi)-y(/li,-y)) (Ai+y) 


y - a 


If we set y = Zq and a = 0, 

y(di,vi)= V(RuZi), 
V(Ai,-y)^ y(7?= VI7,Z = 0), 


(A.7) 


(A.8) 

(A.9) 


and 

Al = 


(/?? + z? - zg) + i ^(R^ + z? - z2)2 + 4Ry^ . (A.IO) 

Thus, from Fqs. A.4 and A.7 (similar to Fq.[^, we obtain a sim¬ 
ple expression for the third integral X at position (di, vi), which 
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is exact in the case of a Stackel potential and explicitly depends 
on the potential, the coordinates, and the velocities. We note that, 
in practice, the intermediate functions h(A) and h(v) do not need 
to be evaluated. 

In summary, the proposed quasi integral in Eq. [^results from 
the exact third integral of the orbits in the Stackel potential of 
Eq. |A.2| which is equal to the true potential on the surfaces A - 
Al and v = vi. The integral is re-evaluated locally for each point 
(T, v) = (/li,vi). 

The case of Stackel potentials with oblate spheroidal coordi¬ 
nates leads to the same equations, but with zo imaginary, Zq < 0 
and y < a. In the zq = limit, a Stackel potential is separable 
in R and z, 

ViR,z)^ViiR) + V 2 iz), 

and 

-Is = [V2(z) - VzCO)] + vhl. 

Thus the vertical energy is an integral of motion. 

In the Zo = 0 limit, Stackel potentials have the form 

ViiO) 

V{R,z)^Vi{r)+^, 

yl 

where - R^ + z^ and 6 is the polar angle 
Table 1), and the integral of motion is 

-zlls ^ Wiie) - Viinm + - l\) . 


(Lynden-Bell|l962 
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